load libraries

library(psych)
library(reticulate)
library(corrplot)
library(lm.beta)
library(boot)
library(plyr)
library(tibble)
library(lme4)
library(lmerTest)
library(NbClust)
library(zoo)
library(ggplot2)
library(ggsci)
library(Hmisc)
library(stringr)


# set virtual environment 
reticulate::use_python("/Users/yluan/.virtualenvs/r_virtualenv/bin/python3.8/", required = TRUE) ## set the path to python under virtual environment 

# load toolbox for generating surrogate maps for spin test
# reticulate::py_install("netneurotools", pip_options="--no-binary='netneurotools'", ignore_installed=TRUE)
# reticulate::py_install("neuromaps", pip_options="--no-binary='neuromaps'", ignore_installed=TRUE)
np  <- import("numpy", convert = F)
nnt <- import("netneurotools")
nntdata <- import("netneurotools")$datasets

nm <- import("neuromaps")
nmnulls <- import("neuromaps")$nulls
nmimgs <- import("neuromaps")$images
nmparc <- import("neuromaps")$parcellate

# convert the schaefer atlas to a gifti file
schaefer = nntdata$fetch_schaefer2018('fslr32k')['200Parcels7Networks']
parc <- nmimgs$dlabel_to_gifti(schaefer)



# lead Schaefer2018_200Parcels_7Networks_order_numeric.txt 
# set the path to schaefer200x7CommunityAffiliation_corrected.txt
sf.assign <- read.table("/Volumes/Local3/Projects/Huashan/scripts/stats/simulation_rmd/schaefer200x7CommunityAffiliation_corrected.txt") 
colors.Schaefer.networks <- c("forestgreen" ,"indianred2" ,"darkgoldenrod1" ,"lemonchiffon1" ,"skyblue3" ,"mediumorchid1" ,"magenta4")[c(3,7,6,5,2,4,1)]

—– PREPARE SIMULATION DATA ——

define numbers of brain ROIs and subjects

n.ROIs <-  200
n.subjects = 150

simulate connectivity matrix

# functional connectivity matrix
set.seed(45) 
fc.mat <- matrix(runif(n.ROIs * n.ROIs, min = -1, max = 1), nrow = n.ROIs, ncol = n.ROIs)
diag(fc.mat) <- 1
fc.mat <- fisherz(fc.mat)
fc.mat[which(fc.mat==Inf,arr.ind = T)] <- NA

# 1000 shuffled functional connectivity matrices
shuffled.fc <- lapply(1:1000, function(x){
  fc.mat.tmp <- matrix(runif(n.ROIs * n.ROIs, min = -1, max = 1), nrow = n.ROIs, ncol = n.ROIs)
  diag(fc.mat.tmp) <- 1
  fc.mat.tmp <- fisherz(fc.mat.tmp)
  fc.mat.tmp[which(fc.mat.tmp==Inf,arr.ind = T)] <- NA
  fc.mat.tmp
})

# Euclidean distance matrix
eu.d <- matrix(runif(n.ROIs * n.ROIs, min = 0, max = 15), nrow = n.ROIs, ncol = n.ROIs)
diag(eu.d) <- NA

simulate subject characteristics

set.seed(45) 

df <- data.frame(SubID = 1:n.subjects,
                 Age.SDM8 = sample(60:90, n.subjects, replace = TRUE),
                 SEX = sample(c("male", "female"), n.subjects, replace = TRUE),
                 DX.AB = sample(c("CN_0", "CN_1","MCI_1", "Dementia_1"), n.subjects, replace = TRUE),
                 pTau181 = runif(n.subjects),
                 AV45.global.mean = runif(n.subjects, min = 0, max = 2.5)) 

df$DX <- mapply(function(x){strsplit(x,"_")[[1]][1]}, df$DX.AB)
df$AV45.bin <- mapply(function(x){as.numeric(strsplit(x,"_")[[1]][2])}, df$DX.AB)

simulate SV2A-PET SUVRs, SV2A-PET w-score and SV2A-PET w-score change rates

# Generate random SV2A PET SUVRs for each ROI for the same subjects
set.seed(45) 

SV2A_PET_SUVRs <- matrix(runif(n.subjects * n.ROIs, min = 0, max = 2), nrow = n.subjects)
colnames(SV2A_PET_SUVRs) <- paste0("SDM8.V", 1:n.ROIs)

# Generate random SV2A PET w-score  for each ROI for the same subjects
SV2A_PET_wscore <- matrix(runif(n.subjects * n.ROIs, min = -2, max = 2), nrow = n.subjects)
colnames(SV2A_PET_wscore) <- paste0("W.SDM8.V", 1:n.ROIs)

# Generate random SV2A PET w-score change rates for each ROI for a subset of subjects
SV2A_PET_wscore_CR <- matrix(runif(25 * n.ROIs, min = -1, max = 0.5), nrow = 25)
colnames(SV2A_PET_wscore_CR) <- paste0("CR.W.SDM8.V", 1:n.ROIs)
SV2A_PET_wscore_CR <- as.data.frame(SV2A_PET_wscore_CR)
SV2A_PET_wscore_CR$SubID <- sample(1:124,25,replace = F)


# Combine the sv2a PET SUVRs and w score with the existing data frame
df <- cbind(df, SV2A_PET_SUVRs, SV2A_PET_wscore)
df <- merge(df, SV2A_PET_wscore_CR, by = "SubID", all.x = T)

# group-average SV2A w-score
df.ad <- subset(df, DX.AB %in% c("CN_1","MCI_1","Dementia_1"))
mean.w.ad <- colMeans(df.ad[,paste0("W.SDM8.V",1:200)])
mean.w.ad.surrogates <- nmnulls$vasa(mean.w.ad,atlas='fsLR',density='32k',
                                            parcellation=parc,seed=as.integer(45))

df.cn0 <- subset(df, DX.AB %in% c("CN_0"))
mean.w.cn0  <- colMeans(df.cn0[,paste0("W.SDM8.V",1:200)])

—– ANALYSES ——

Analysis 1: the association between sv2a covariance and functional connectivity matrix

# generate the sv2s covariance matrix in abeta positive group
w.cov.ad <- cor(as.matrix(df.ad[,paste0("W.SDM8.V",1:200)]),method = "spearman")
w.cov.z.ad <- fisherz(w.cov.ad)
w.cov.z.ad[which(w.cov.z.ad==Inf)] <- NA
corrplot(w.cov.z.ad, diag = FALSE, tl.pos = "n", tl.cex = 0.5, method = "color",is.corr = F)

# generate the sv2s covariance matrix in cognitively normal abeta negative group
w.cov.cn0 <- cor(as.matrix(df.cn0[,paste0("W.SDM8.V",1:200)]),method = "spearman")
w.cov.z.cn0 <- fisherz(w.cov.cn0)
w.cov.z.cn0[which(w.cov.z.cn0==Inf)] <- NA
corrplot(w.cov.z.cn0, diag = FALSE, tl.pos = "n", tl.cex = 0.5, method = "color",is.corr = F)

# prepare data fram for analyses
w.cov.vecterized.ad <- w.cov.z.ad[lower.tri(w.cov.z.ad)]
w.cov.vecterized.cn0 <- w.cov.z.cn0[lower.tri(w.cov.z.cn0)]
fc.verterized <- fc.mat[lower.tri(fc.mat)]
df.stat.cov <- data.frame(w.cov.ad= w.cov.vecterized.ad,
                          w.cov.cn0=w.cov.vecterized.cn0,
                          fc=fc.verterized)

## test for Abeta+ group
mod.lm.ad <- lm(data = df.stat.cov,w.cov.ad~fc)
mod.lm.ad <- summary(lm.beta(mod.lm.ad))
beta.val.ad <- mod.lm.ad$coefficients["fc","Standardized"]
mod.lm.ad
## 
## Call:
## lm(formula = w.cov.ad ~ fc, data = df.stat.cov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.46136 -0.06379  0.00019  0.06460  0.41015 
## 
## Coefficients:
##               Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept) -3.695e-04           NA  6.820e-04  -0.542    0.588
## fc           7.572e-06    7.195e-05  7.461e-04   0.010    0.992
## 
## Residual standard error: 0.0962 on 19898 degrees of freedom
## Multiple R-squared:  5.176e-09,  Adjusted R-squared:  -5.025e-05 
## F-statistic: 0.000103 on 1 and 19898 DF,  p-value: 0.9919
beta.shuffled.ad <- mapply(function(x){
  fc.mat.shuffled <- shuffled.fc[[x]]
  df.cov.shuffled <- data.frame(w.cov.ad= w.cov.vecterized.ad,
                                w.cov.cn0=w.cov.vecterized.cn0,
                                fc=fc.mat.shuffled[lower.tri(fc.mat.shuffled)])
  mod.lm.shuffled <- lm(data = df.cov.shuffled,w.cov.ad~fc)
  mod.lm.shuffled  <- summary(lm.beta(mod.lm.shuffled))
  beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
}, 1:1000)
p.val.shuffled.ad <- length(which(beta.shuffled.ad > beta.val.ad))/length(beta.shuffled.ad)

## test for CN abeta- group
mod.lm.cn0 <- lm(data = df.stat.cov,w.cov.cn0~fc)
mod.lm.cn0 <- summary(lm.beta(mod.lm.cn0))
beta.val.cn0 <- mod.lm.cn0$coefficients["fc","Standardized"]
mod.lm.cn0
## 
## Call:
## lm(formula = w.cov.cn0 ~ fc, data = df.stat.cov)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.71767 -0.11431 -0.00083  0.11342  0.72274 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)  0.001930           NA   0.001221   1.580    0.114
## fc          -0.001896    -0.010061   0.001336  -1.419    0.156
## 
## Residual standard error: 0.1723 on 19898 degrees of freedom
## Multiple R-squared:  0.0001012,  Adjusted R-squared:  5.098e-05 
## F-statistic: 2.014 on 1 and 19898 DF,  p-value: 0.1558
beta.shuffled.cn0 <- mapply(function(x){
  fc.mat.shuffled <- shuffled.fc[[x]]
  df.cov.shuffled <- data.frame(w.cov.ad= w.cov.vecterized.ad,
                                w.cov.cn0=w.cov.vecterized.cn0,
                                fc=fc.mat.shuffled[lower.tri(fc.mat.shuffled)])
  mod.lm.shuffled <- lm(data = df.cov.shuffled,w.cov.cn0~fc)
  mod.lm.shuffled  <- summary(lm.beta(mod.lm.shuffled))
  beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
}, 1:1000)
p.val.shuffled.cn0 <- length(which(beta.shuffled.cn0 > beta.val.cn0))/length(beta.shuffled.cn0)


## comparison between Abeta+ group and CN abeta- group with bootstrapping
### define function
sv2a.cov.z.boot <- function(mat,ind){
  sv2a.mat <- as.matrix(mat[ind,paste0("W.SDM8.V",1:200)])
  w.cov.boot <- cor(sv2a.mat,method = "spearman")
  w.cov.z.boot <- fisherz(w.cov.boot)
  w.cov.z.boot[which(w.cov.z.boot==Inf)] <- NA
  
  w.cov.vecterized.boot<- w.cov.z.boot[lower.tri(w.cov.z.boot)]
  fc.verterized <- fc.mat[lower.tri(fc.mat)]
  
  df.cov.boot <- data.frame(w.cov.boot= w.cov.vecterized.boot,
                            fc=fc.verterized)

  mod.lm <- lm(data = df.cov.boot,w.cov.boot~fc)
  mod.lm <- summary(lm.beta(mod.lm))
  beta.val <- mod.lm$coefficients["fc","Standardized"]

  return(beta.val)
}


beta.ad.boot <- boot(data=df.ad,statistic=sv2a.cov.z.boot,R=1000)
beta.cn0.boot <- boot(data=df.cn0,statistic=sv2a.cov.z.boot,R=1000)

mod.t <- t.test(beta.ad.boot$t, beta.cn0.boot$t)
mod.t
## 
##  Welch Two Sample t-test
## 
## data:  beta.ad.boot$t and beta.cn0.boot$t
## t = 31.072, df = 1998, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.006763386 0.007674673
## sample estimates:
##     mean of x     mean of y 
## -9.940842e-06 -7.228971e-03

Analysis 2: network spreading Model (Shafiei, et al, Biol Psychiatry 2020)

# define functions 
get.network.spreading.beta.spin <- function(mean.sv2a.w,surrogate.maps,fc.thre){
  mapply(function(x){
    mean.w.ad.surrog.curren <- surrogate.maps[,x]
    
    fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
    fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]

    fc.w.sv2a.ad <- mapply(function(k){
      fc.tmp <- fc.mat
      fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
      fc.w.tmp <- mean(mean.w.ad.surrog.curren[-k]*fc.tmp[-k,k],na.rm=T)
      fc.w.tmp
    }, 1:200)
    
    df.stat.spin <- data.frame(fc.wei.w.ad= fc.w.sv2a.ad,
                               mean.w = mean.sv2a.w)
    
    mod.lm <- lm(data = df.stat.spin,mean.w~fc.wei.w.ad)
    mod.lm <- summary(lm.beta(mod.lm))
    beta.val.spin <- mod.lm$coefficients[2,2]
    
    return(beta.val.spin)
  }, 1:1000)
}

get.network.spreading.beta.shuffled <- function(mean.sv2a.w, fc.thre){
  mapply(function(x){
    fc.mat.shuffled <- as.matrix(shuffled.fc[[x]])
    
    fc.shuffled.vals <- fc.mat.shuffled[as.vector(lower.tri(fc.mat.shuffled))]
    fc.shuffled.thre.val <- fc.shuffled.vals[order(fc.shuffled.vals)][length(fc.shuffled.vals)*fc.thre+1]
    
    fc.w.sv2a.ad.shuffled <- mapply(function(k){
      fc.tmp <- fc.mat.shuffled
      fc.tmp[which(fc.tmp<fc.shuffled.thre.val)] <- NA
      fc.w.tmp <- mean(mean.sv2a.w[-k]*fc.tmp[-k,k],na.rm=T)
      fc.w.tmp
    }, 1:200)
    
    df.stat.shuffled <- data.frame(fc.wei.w.ad= fc.w.sv2a.ad.shuffled,
                                   mean.w = mean.sv2a.w)
    
    mod.lm.shuffled <- lm(data = df.stat.shuffled,mean.w~fc.wei.w.ad)
    mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
    beta.val.shuffled <- mod.lm.shuffled$coefficients[2,2]
 
    return(beta.val.shuffled)
  }, 1:1000)
}



# analyses 
fc.thre.all <- seq(0,1,0.25)[1:4] # set different functional network sparsity

beta.null.all <- list()
beta.empirical.all <- c()
p.spin.all <- c()
p.shuffled.all <- c()

for (i in 1:4) {
  ## empirical model
  fc.thre <- fc.thre.all[i]
  fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
  fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]
  
  
  fc.w.sv2a.tmp <- mapply(function(x){
    fc.tmp <- fc.mat
    fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
    fc.w.tmp <- mean(mean.w.ad[-x]*fc.tmp[-x,x],na.rm=T)
    c(fc.w.tmp)
  }, 1:200)
  
  df.stat.tmp <- data.frame(fc.wei.w = fc.w.sv2a.tmp,mean.w = mean.w.ad)
  mod.lm <- lm(data = df.stat.tmp,mean.w~fc.wei.w)
  mod.lm <- summary(lm.beta(mod.lm))
  mod.lm
  beta.val.empirical <- mod.lm$coefficients["fc.wei.w","Standardized"]
  p.val.empirical <- mod.lm$coefficients["fc.wei.w","Pr(>|t|)"]
  
  ## generate null distribution of beta values
  null.beta.spin <- get.network.spreading.beta.spin(mean.w.ad,mean.w.ad.surrogates, fc.thre)
  null.beta.shuffled <- get.network.spreading.beta.shuffled(mean.w.ad, fc.thre)
  
  df.null.beta <- data.frame(null = c(rep("spin",1000),rep("rewired",1000)),
                             beta = c(null.beta.spin,null.beta.shuffled))
  p.val.spin <- length(which(null.beta.spin > beta.val.empirical))/length(null.beta.spin)
  p.val.shuffled <- length(which(null.beta.shuffled > beta.val.empirical))/length(null.beta.shuffled)

  beta.null.all[[i]] <- df.null.beta
  beta.empirical.all[i] <- beta.val.empirical
  p.spin.all[i] <- p.val.spin
  p.shuffled.all[i] <- p.val.shuffled
}

beta.empirical.all <- data.frame(beta.empirical.all)
beta.empirical.all$x.point <- c(1,2,3,4)


df.beta.null.all <- data.frame(Sparsity=c(rep("100%",2000),rep("75%",2000),rep("50%",2000), rep("25%",2000)),
                               null.model=rep(c(rep("spin",1000),rep("shuffled",1000)),4),
                               beta.val= c(beta.null.all[[1]][,2],beta.null.all[[2]][,2],beta.null.all[[3]][,2],beta.null.all[[4]][,2]))
df.beta.null.all$Sparsity <- factor(df.beta.null.all$Sparsity, levels = c("100%","75%","50%","25%"))

y.p.val <- min(df.beta.null.all$beta.val)

significance <- c(p.spin.all[1],p.shuffled.all[1],p.spin.all[2],p.shuffled.all[2],p.spin.all[3],p.shuffled.all[3],p.spin.all[4],p.shuffled.all[4])
significance <- as.data.frame(significance)
significance$null.model <- c(rep(c("spin","shuffled"),4))
significance$Sparsity <- c(rep(c("100%","75%","50%","25%"),each = 2))
significance$y.postion <- y.p.val

ggplot(data=df.beta.null.all)+
  aes(x=Sparsity, y=beta.val)+
  geom_point(position = position_jitterdodge(0.65),shape=16, size=0.35, alpha=0.55,aes(fill = factor(null.model),color = factor(null.model)))+
  geom_boxplot(notch = F,color="black", alpha=0,aes(fill = factor(null.model),color = factor(null.model)),linewidth=0.25)+
  ylab("Regression-derived β-value")+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color= " black", fill = NA),
        axis.text.x = element_text(size=12),
        axis.text.y = element_text(size=10),
        axis.title.y = element_text(size=14))+
  scale_color_manual(values=c("#79AF97FF", "#6A6599FF"))+
  geom_point(data=beta.empirical.all,aes(x=x.point,y=beta.empirical.all),colour="#F56F5C",size=2)+
  geom_text(data=significance,aes(label = ifelse(significance<0.001, "***",ifelse(significance<0.01,"**",ifelse(significance<0.05,"*",""))), 
                                  group = null.model,y = y.postion-0.05), 
            position = position_dodge(width = .75), vjust = 0.5,size = 16 / .pt)

# subject-level analyses
for (i in 1:4) {
  fc.thre <- fc.thre.all[i]
  print(paste0("The functional network sparsity at ", fc.thre))
  fc.vals <- fc.mat[as.vector(lower.tri(fc.mat))]
  fc.thre.val <- fc.vals[order(fc.vals)][length(fc.vals)*fc.thre+1]
  
  beta.network.spreading.sub <- ddply(df.ad,.(SubID),function(x){
    w.ad.sub <- x[,paste0("W.SDM8.V",1:200)]
    w.ad.sub <- as.vector(as.matrix(w.ad.sub))
    
    fc.w.sv2a.sub <- mapply(function(x){
      fc.tmp <- fc.mat
      fc.tmp[which(fc.tmp<fc.thre.val)] <- NA
      fc.w.tmp <- mean(w.ad.sub[-x]*fc.tmp[-x,x],na.rm=T)
      fc.w.tmp
    }, 1:200)

    df.stat.ad.sub <- data.frame(fc.wei.w.sub = fc.w.sv2a.sub,mean.w = w.ad.sub)
    
    mod.lm <- lm(data = df.stat.ad.sub,mean.w~fc.wei.w.sub)
    mod.lm <- summary(lm.beta(mod.lm))
    beta.val.sub <- mod.lm$coefficients[2,2]
    beta.val.sub
  })
  
  beta.network.spreading.sub$variable <- "beta"
  
  t.mod <- t.test(beta.network.spreading.sub$V1)
  print(t.mod)
}
## [1] "The functional network sparsity at 0"
## 
##  One Sample t-test
## 
## data:  beta.network.spreading.sub$V1
## t = 0.71183, df = 112, p-value = 0.478
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.00832983  0.01767095
## sample estimates:
##   mean of x 
## 0.004670561 
## 
## [1] "The functional network sparsity at 0.25"
## 
##  One Sample t-test
## 
## data:  beta.network.spreading.sub$V1
## t = -4.3518, df = 112, p-value = 2.996e-05
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.03922860 -0.01468265
## sample estimates:
##   mean of x 
## -0.02695563 
## 
## [1] "The functional network sparsity at 0.5"
## 
##  One Sample t-test
## 
## data:  beta.network.spreading.sub$V1
## t = -6.8285, df = 112, p-value = 4.66e-10
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.05591255 -0.03076263
## sample estimates:
##   mean of x 
## -0.04333759 
## 
## [1] "The functional network sparsity at 0.75"
## 
##  One Sample t-test
## 
## data:  beta.network.spreading.sub$V1
## t = -6.3251, df = 112, p-value = 5.346e-09
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.05030069 -0.02630383
## sample estimates:
##   mean of x 
## -0.03830226

Analysis 3: epicenter connectivity-based sv2a propagation in pooled abeta+ subjects

# define functions
get.epicenter.by.rank <- function(mean.sv2a.w){
  rank.mean.w <- as.vector(rank(mean.sv2a.w))
  fc.wei.w <- mapply(function(k){
    fc.w.tmp <- mean(mean.sv2a.w[-k]*fc.mat[-k,k])
    fc.w.tmp
  }, 1:200)
  rank.fc.wei.w <- rank(fc.wei.w)
  mean.rank.current <- (rank.mean.w+rank.fc.wei.w)/2
  
  ## generated null distribution of mean ranks
  surrogate.map <- nmnulls$vasa(mean.sv2a.w,atlas='fsLR',density='32k',parcellation=parc,seed=as.integer(45))
  
  mean.rank.null <- mapply(function(j){
    mean.w.surrog <- surrogate.map[,j]
    rank.mean.w.surrog <- as.vector(rank(mean.w.surrog))
    
    fc.wei.w.surrog<- mapply(function(k){
      fc.w.tmp <- mean(mean.w.surrog[-k]*fc.mat[-k,k])
      fc.w.tmp
    }, 1:200)
    rank.fc.wei.w.surrog <- rank(fc.wei.w.surrog)
    mean.rank <- (rank.mean.w.surrog+rank.fc.wei.w.surrog)/2
    
  }, 1:1000)
  
  rank.p.vals <- mapply(function(roi){
    p.roi <- length(which(mean.rank.current[roi] > mean.rank.null[roi,]))/length(mean.rank.null[roi,])
  }, 1:200)
  
  epicenter.roi <- rep(0,n.ROIs)
  epicenter.roi[which(rank.p.vals<0.05)] <- 1
  
  epicenter.roi
}


get.epi.pred.beta.shuffled <- function(epicenter.rois, mean.sv2a.w){
  mapply(function(j){
    fc.mat.shuffled <- shuffled.fc[[j]]
    
    epi.fc.mat.shuffled <- as.matrix(fc.mat.shuffled[epicenter.rois,])
    epi.fc.mat.shuffled[which(epi.fc.mat.shuffled==Inf)] <- NA
    mean.epi.fc.mat.shuffled <- colMeans(epi.fc.mat.shuffled,na.rm = T)
    
    df.test.ad.shuffled <- data.frame(epi.fc = mean.epi.fc.mat.shuffled,
                                      sv2a.wscore=mean.sv2a.w)
    mod.lm.shuffled <- lm(data = df.test.ad.shuffled,sv2a.wscore~epi.fc)
    mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
    beta.val <- mod.lm.shuffled$coefficients[2,2]
    beta.val
  }, 1:1000)

}



# Epicenter connectivity-based prediciton of group-level sv2a w score in pooled abeta+ subjects
## get group-level epicenter ROIs
epicenter.ad <- get.epicenter.by.rank(mean.w.ad)
epicenter.ad <- which(epicenter.ad==1)

## epicenter-based prediction
epi.fc.mat.ad <- fc.mat[epicenter.ad,]
mean.epi.fc.ad <- colMeans(epi.fc.mat.ad,na.rm = T)

## prediction
beta.shuffled.ad <- get.epi.pred.beta.shuffled(epicenter.ad,mean.w.ad)

df.stat.epicenter.fc.prediction <- data.frame(epicenter.fc = mean.epi.fc.ad, mean.w = mean.w.ad)
mod.lm <- lm(data = df.stat.epicenter.fc.prediction, mean.w~epicenter.fc)
mod.lm <- summary(lm.beta(mod.lm))
mod.lm
## 
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.epicenter.fc.prediction)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.35163 -0.09284  0.00231  0.07666  0.31812 
## 
## Coefficients:
##                Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)  -0.0001211           NA  0.0082469  -0.015    0.988
## epicenter.fc -0.0250385   -0.0617498  0.0287614  -0.871    0.385
## 
## Residual standard error: 0.1166 on 198 degrees of freedom
## Multiple R-squared:  0.003813,   Adjusted R-squared:  -0.001218 
## F-statistic: 0.7579 on 1 and 198 DF,  p-value: 0.385
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]

p.val.shuffled <- length(which(beta.shuffled.ad < beta.val))/length(beta.shuffled.ad)


## plot
x.range <- range(df.stat.epicenter.fc.prediction$epicenter.fc)
y.range <- range(df.stat.epicenter.fc.prediction$mean.w)
tibble(Predicted = df.stat.epicenter.fc.prediction$epicenter.fc, Response = df.stat.epicenter.fc.prediction$mean.w, nets = as.factor(sf.assign$V1)) %>%
  ggplot(aes(x=Predicted, y = Response))+
  geom_point(shape=21, size=2,aes(fill=factor(nets)))+
  geom_smooth(method = "lm",color = "grey55")+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color= " black", fill = NA),
        aspect.ratio = 1)+
  xlab("Functional connectivity to epicenter")+
  ylab("SV2A-PET W-score")+
  scale_fill_manual(values=colors.Schaefer.networks,labels=c("Visual","Motor","DAN","VAN","Limbic","FPCN","DMN"),name="")+
  scale_y_continuous(labels = scales::number_format(accuracy = 0.01))+
  annotate("text", x =  min(x.range)+2/5*diff(x.range), y = max(y.range)+1/8*diff(y.range),
           label = paste0(paste0("italic(β)==", round(beta.val,3)), "*\',\'",
                          "~", expression(italic(p)[italic(rewired)]),ifelse(p.val.shuffled>0.001,paste0("==",round(p.val.shuffled,3)),"<0.001")),
           color = "black",parse=T, size=5)
## `geom_smooth()` using formula = 'y ~ x'

# Epicenter connectivity-based prediciton of subject-level sv2a w score in pooled abeta+ subjects
## define subject-specific epicenter
 df.epicenters.subject <- ddply(df.ad,.(SubID), function(x){
      mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
      epicenters.sub <- get.epicenter.by.rank(mean.w.sub)
      epicenters.sub
    })
colnames(df.epicenters.subject)[2:201] <- paste0("Epicenter.bin.V",1:200)
df.ad.epicenter <- merge(df.ad,df.epicenters.subject, by ="SubID")

## subject-level prediction
df.stat.sub <- ddply(df.ad.epicenter,.(SubID), function(x){

  mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
  epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
  
  epi.fc.mat.sub <- fc.mat[epicenters.sub,]
  mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
  
  fc.quantile.thre <-  quantile(mean.epi.fc.sub)
  q1.rois <- which(mean.epi.fc.sub<fc.quantile.thre[2])
  q2.rois <- which(mean.epi.fc.sub<fc.quantile.thre[3]&mean.epi.fc.sub>=fc.quantile.thre[2])
  q3.rois <- which(mean.epi.fc.sub<fc.quantile.thre[4]&mean.epi.fc.sub>=fc.quantile.thre[3])
  q4.rois <- which(mean.epi.fc.sub>=fc.quantile.thre[4])
  
  fc.q1 <- mean(mean.epi.fc.sub[q1.rois])
  fc.q2 <- mean(mean.epi.fc.sub[q2.rois])
  fc.q3 <- mean(mean.epi.fc.sub[q3.rois])
  fc.q4 <- mean(mean.epi.fc.sub[q4.rois])

  w.q1 <- mean(mean.w.sub[q1.rois])
  w.q2 <- mean(mean.w.sub[q2.rois])
  w.q3 <- mean(mean.w.sub[q3.rois])
  w.q4 <- mean(mean.w.sub[q4.rois])
  
  epi.eu.mat.sub <- eu.d[epicenters.sub,]
  mean.epi.eu.d.sub <- colMeans(epi.eu.mat.sub,na.rm = T)
  eu.q1 <- mean(mean.epi.eu.d.sub[q1.rois])
  eu.q2 <- mean(mean.epi.eu.d.sub[q2.rois])
  eu.q3 <- mean(mean.epi.eu.d.sub[q3.rois])
  eu.q4 <- mean(mean.epi.eu.d.sub[q4.rois])
  
  df.fc.sub <- data.frame(Quantile = c("Q1","Q2","Q3","Q4"),
                          Quantile.num = 1:4,
                          FC.quantile = c(fc.q1,fc.q2,fc.q3,fc.q4),
                          mean.sv2a = c(w.q1,w.q2,w.q3,w.q4),
                          EU.dist = c(eu.q1,eu.q2,eu.q3,eu.q4))
  
})

df.stat.sub <- merge(df.ad.epicenter,df.stat.sub, by ="SubID")
mod.lmer <- lmer(data = df.stat.sub,mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB +EU.dist +(1|SubID)) ;summary(mod.lmer)
## boundary (singular) fit: see help('isSingular')
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +  
##     (1 | SubID)
##    Data: df.stat.sub
## 
## REML criterion at convergence: -329
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0007 -0.6760 -0.0116  0.6323  3.0220 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.00000  0.0000  
##  Residual             0.02575  0.1605  
## Number of obs: 452, groups:  SubID, 113
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)     -3.515e-01  2.815e-01  4.450e+02  -1.249   0.2124  
## FC.quantile      2.079e-02  2.761e-02  4.450e+02   0.753   0.4520  
## SEXmale          2.864e-02  1.562e-02  4.450e+02   1.833   0.0674 .
## Age.SDM8        -5.009e-04  8.678e-04  4.450e+02  -0.577   0.5641  
## DX.ABDementia_1 -1.037e-03  1.799e-02  4.450e+02  -0.058   0.9541  
## DX.ABMCI_1      -1.840e-02  1.917e-02  4.450e+02  -0.960   0.3375  
## EU.dist          5.086e-02  3.761e-02  4.450e+02   1.352   0.1770  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.057                                   
## SEXmale     -0.002 -0.010                            
## Age.SDM8    -0.121 -0.020  0.072                     
## DX.ABDmnt_1  0.062  0.003  0.041 -0.067              
## DX.ABMCI_1   0.040 -0.003 -0.187  0.053  0.428       
## EU.dist     -0.972  0.062 -0.038 -0.111 -0.077 -0.076
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## plot
y.range <- range(df.stat.sub$mean.sv2a)
ggplot(data= df.stat.sub)+
  aes(x=Quantile, y = mean.sv2a)+
  geom_point( size=2, color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
  geom_line(color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
  geom_smooth(method = "lm",color = "grey55",aes(x=Quantile.num, y = mean.sv2a))+
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color= " black", fill = NA),
        axis.text.x = element_text(size=10),
        axis.text.y = element_text(size=10))+
  xlab(str_wrap("Functional connectivity quantile to subject-specific epicenter", width = 40))+
  ylab(str_wrap("SV2A-PET W-score", width = 30))
## `geom_smooth()` using formula = 'y ~ x'

Analysis 4: epicenter connectivity-based sv2a propagation in subtypes of abeta+ subjects

# hierachical clustering to determine three sv2a w-score subtypes
df.cluster.ad <- df.ad.epicenter
df.cluster.ad.w <- df.cluster.ad[,paste0("W.SDM8.V",1:200)]
df.cluster.ad.w <- scale(t(df.cluster.ad.w))
df.cluster.ad.w <- t(df.cluster.ad.w)

d <- dist(df.cluster.ad.w,method = 'euclidean')
fit.hc <- hclust(d, method = "ward.D")
clusters <- cutree(fit.hc,k=3)
df.ad.epicenter$clusters <- clusters


# group-level epicenter connectivity-based prediction of sv2a w-score for each subtypese
for (i in 1:3) {
  print(paste0("Analysis for cluster ", i))
  df.ad.epicenter.cluster <- subset(df.ad.epicenter, clusters==i)
  mean.w.ad.cluster <- colMeans(df.ad.epicenter.cluster[,paste0("W.SDM8.V",1:200)])
  epicenter.ad.cluster <- get.epicenter.by.rank(mean.w.ad.cluster)
  epicenter.ad.cluster <- which(epicenter.ad.cluster==1)
  
  
  ## epicenter-based prediction
  epi.fc.mat.ad.cluster <- fc.mat[epicenter.ad.cluster,]
  mean.epi.fc.ad.cluster <- colMeans(epi.fc.mat.ad.cluster,na.rm = T)
  
  ## prediction
  beta.shuffled.ad.cluster <- get.epi.pred.beta.shuffled(epicenter.ad.cluster, mean.w.ad.cluster)
  
  df.stat.cluster.epicenter.fc.prediction <- data.frame(epicenter.fc = mean.epi.fc.ad.cluster, mean.w = mean.w.ad.cluster)
  mod.lm <- lm(data = df.stat.cluster.epicenter.fc.prediction, mean.w~epicenter.fc)
  mod.lm <- summary(lm.beta(mod.lm))
  print(mod.lm)
  beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
  p.val.shuffled <- length(which(beta.shuffled.ad.cluster < beta.val))/length(beta.shuffled.ad.cluster)
  
  
  ## plot
  x.range <- range(df.stat.cluster.epicenter.fc.prediction$epicenter.fc)
  y.range <- range(df.stat.cluster.epicenter.fc.prediction$mean.w)
  tibble(Predicted = df.stat.cluster.epicenter.fc.prediction$epicenter.fc, Response = df.stat.cluster.epicenter.fc.prediction$mean.w, nets = as.factor(sf.assign$V1)) %>%
    ggplot(aes(x=Predicted, y = Response))+
    geom_point(shape=21, size=2,aes(fill=factor(nets)))+
    geom_smooth(method = "lm",color = "grey55")+
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          panel.border = element_rect(color= " black", fill = NA),
          aspect.ratio = 1)+
    xlab("Functional connectivity to epicenter")+
    ylab("SV2A-PET W-score")+
    ggtitle(paste0("Epicenter connectivity-based prediction of SV2A for subtype ",i))+
    scale_fill_manual(values=colors.Schaefer.networks,labels=c("Visual","Motor","DAN","VAN","Limbic","FPCN","DMN"),name="")+
    scale_y_continuous(labels = scales::number_format(accuracy = 0.01))+
    annotate("text", x =  min(x.range)+2/5*diff(x.range), y = max(y.range)+1/8*diff(y.range),
             label = paste0(paste0("italic(β)==", round(beta.val,3)), "*\',\'",
                            "~", expression(italic(p)[italic(rewired)]),ifelse(p.val.shuffled>0.001,paste0("==",round(p.val.shuffled,3)),"<0.001")),
             color = "black",parse=T, size=5)
  
  ## sliding window analyses
  df.ad.epicenter.cluster$mean.sv2a <- rowMeans(df.ad.epicenter.cluster[,paste0("W.SDM8.V",1:200)])
  sub.order.cluster <- order(df.ad.epicenter.cluster$mean.sv2a,decreasing = T)
  df.ad.epicenter.cluster.slide <- df.ad.epicenter.cluster[sub.order.cluster,]
  num.sub <- length(sub.order.cluster)
  win.wid <- 10 # set the width of windows at 10 subjects
  win.num <- 9  # set the number of windows at 10 windows
  sep <- win.wid-((win.wid*win.num-num.sub)/(win.num-1))
  windows.assign <- rollapply(1:num.sub, win.wid, by = sep, c)

  
  for (i.win in 1:win.num) {
    df.current.window <- df.ad.epicenter.cluster.slide[as.vector(windows.assign[i.win,]),]
    mean.w.current.window <- colMeans(df.current.window[,paste0("W.SDM8.V",1:200)])
    epicenters.current.window <- epicenter.ad.cluster
    epi.fc.current.window <- fc.mat[epicenters.current.window,]
    
    df.test.current.window <- data.frame(epicenter.fc = colMeans(epi.fc.current.window,na.rm = T),
                                         mean.w = mean.w.current.window)
    
    beta.shuffled.current.window <- get.epi.pred.beta.shuffled(epicenters.current.window,mean.w.current.window)
    mod.lm <- lm(data = df.test.current.window,mean.w~epicenter.fc)
    mod.lm <- summary(lm.beta(mod.lm))
    beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
    p.val.shuffled <- length(which(beta.shuffled.current.window< beta.val))/length(beta.shuffled.current.window)
    r2 <- mod.lm$adj.r.squared
    print(paste0("Sliding window analyses for cluster ", i , " window ", i.win))
    print(paste0("β=", round(beta.val,3) , ", p.rewired=", p.val.shuffled,", r2=",round(r2,3)))
  }
}
## [1] "Analysis for cluster 1"
## 
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.60106 -0.11353  0.00104  0.14385  0.42430 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)  0.007844           NA   0.014211   0.552    0.582
## epicenter.fc 0.091785     0.115261   0.056215   1.633    0.104
## 
## Residual standard error: 0.201 on 198 degrees of freedom
## Multiple R-squared:  0.01329,    Adjusted R-squared:  0.008302 
## F-statistic: 2.666 on 1 and 198 DF,  p-value: 0.1041
## 
## [1] "Sliding window analyses for cluster 1 window 1"
## [1] "β=0.051, p.rewired=0.73, r2=-0.002"
## [1] "Sliding window analyses for cluster 1 window 2"
## [1] "β=0.071, p.rewired=0.828, r2=0"
## [1] "Sliding window analyses for cluster 1 window 3"
## [1] "β=0.036, p.rewired=0.698, r2=-0.004"
## [1] "Sliding window analyses for cluster 1 window 4"
## [1] "β=0.073, p.rewired=0.849, r2=0"
## [1] "Sliding window analyses for cluster 1 window 5"
## [1] "β=0.055, p.rewired=0.776, r2=-0.002"
## [1] "Sliding window analyses for cluster 1 window 6"
## [1] "β=0.102, p.rewired=0.933, r2=0.005"
## [1] "Sliding window analyses for cluster 1 window 7"
## [1] "β=0.125, p.rewired=0.959, r2=0.011"
## [1] "Sliding window analyses for cluster 1 window 8"
## [1] "β=0.096, p.rewired=0.912, r2=0.004"
## [1] "Sliding window analyses for cluster 1 window 9"
## [1] "β=0.023, p.rewired=0.619, r2=-0.005"
## [1] "Analysis for cluster 2"
## 
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -1.06040 -0.24926 -0.02737  0.23933  0.94678 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)  0.004828           NA   0.025959   0.186    0.853
## epicenter.fc 0.116755     0.077688   0.106482   1.096    0.274
## 
## Residual standard error: 0.3662 on 198 degrees of freedom
## Multiple R-squared:  0.006035,   Adjusted R-squared:  0.001015 
## F-statistic: 1.202 on 1 and 198 DF,  p-value: 0.2742
## 
## [1] "Sliding window analyses for cluster 2 window 1"
## [1] "β=0.055, p.rewired=0.785, r2=-0.002"
## [1] "Sliding window analyses for cluster 2 window 2"
## [1] "β=0.068, p.rewired=0.831, r2=0"
## [1] "Sliding window analyses for cluster 2 window 3"
## [1] "β=0.047, p.rewired=0.733, r2=-0.003"
## [1] "Sliding window analyses for cluster 2 window 4"
## [1] "β=0.07, p.rewired=0.817, r2=0"
## [1] "Sliding window analyses for cluster 2 window 5"
## [1] "β=0.067, p.rewired=0.827, r2=-0.001"
## [1] "Sliding window analyses for cluster 2 window 6"
## [1] "β=0.024, p.rewired=0.603, r2=-0.004"
## [1] "Sliding window analyses for cluster 2 window 7"
## [1] "β=0.052, p.rewired=0.748, r2=-0.002"
## [1] "Sliding window analyses for cluster 2 window 8"
## [1] "β=0.042, p.rewired=0.687, r2=-0.003"
## [1] "Sliding window analyses for cluster 2 window 9"
## [1] "β=0.07, p.rewired=0.829, r2=0"
## [1] "Analysis for cluster 3"
## 
## Call:
## lm(formula = mean.w ~ epicenter.fc, data = df.stat.cluster.epicenter.fc.prediction)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -0.7511 -0.1424 -0.0174  0.1561  0.5649 
## 
## Coefficients:
##              Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)  -0.01089           NA    0.01549  -0.703    0.483
## epicenter.fc  0.00197      0.00277    0.05054   0.039    0.969
## 
## Residual standard error: 0.2189 on 198 degrees of freedom
## Multiple R-squared:  7.673e-06,  Adjusted R-squared:  -0.005043 
## F-statistic: 0.001519 on 1 and 198 DF,  p-value: 0.9689
## 
## [1] "Sliding window analyses for cluster 3 window 1"
## [1] "β=-0.054, p.rewired=0.211, r2=-0.002"
## [1] "Sliding window analyses for cluster 3 window 2"
## [1] "β=-0.09, p.rewired=0.103, r2=0.003"
## [1] "Sliding window analyses for cluster 3 window 3"
## [1] "β=-0.076, p.rewired=0.144, r2=0.001"
## [1] "Sliding window analyses for cluster 3 window 4"
## [1] "β=0.063, p.rewired=0.826, r2=-0.001"
## [1] "Sliding window analyses for cluster 3 window 5"
## [1] "β=0.124, p.rewired=0.965, r2=0.01"
## [1] "Sliding window analyses for cluster 3 window 6"
## [1] "β=0.118, p.rewired=0.951, r2=0.009"
## [1] "Sliding window analyses for cluster 3 window 7"
## [1] "β=0.057, p.rewired=0.778, r2=-0.002"
## [1] "Sliding window analyses for cluster 3 window 8"
## [1] "β=-0.071, p.rewired=0.162, r2=0"
## [1] "Sliding window analyses for cluster 3 window 9"
## [1] "β=-0.063, p.rewired=0.175, r2=-0.001"
# subject-level epicenter connectivity-based prediction of sv2a w-score for each subtypese
for (i in 1:3) {
  print(paste0("Subject-level analysis for cluster ", i))
  df.ad.epicenter.cluster <- subset(df.ad.epicenter, clusters==i)
  
  df.stat.sub.cluster <- ddply(df.ad.epicenter.cluster,.(SubID), function(x){
    
    mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
    epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
    
    epi.fc.mat.sub <- fc.mat[epicenters.sub,]
    mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
    
    fc.quantile.thre <-  quantile(mean.epi.fc.sub)
    q1.rois <- which(mean.epi.fc.sub<fc.quantile.thre[2])
    q2.rois <- which(mean.epi.fc.sub<fc.quantile.thre[3]&mean.epi.fc.sub>=fc.quantile.thre[2])
    q3.rois <- which(mean.epi.fc.sub<fc.quantile.thre[4]&mean.epi.fc.sub>=fc.quantile.thre[3])
    q4.rois <- which(mean.epi.fc.sub>=fc.quantile.thre[4])
    
    fc.q1 <- mean(mean.epi.fc.sub[q1.rois])
    fc.q2 <- mean(mean.epi.fc.sub[q2.rois])
    fc.q3 <- mean(mean.epi.fc.sub[q3.rois])
    fc.q4 <- mean(mean.epi.fc.sub[q4.rois])
    
    w.q1 <- mean(mean.w.sub[q1.rois])
    w.q2 <- mean(mean.w.sub[q2.rois])
    w.q3 <- mean(mean.w.sub[q3.rois])
    w.q4 <- mean(mean.w.sub[q4.rois])
    
    
    epi.eu.mat.sub <- eu.d[epicenters.sub,]
    mean.epi.eu.d.sub <- colMeans(epi.eu.mat.sub,na.rm = T)
    eu.q1 <- mean(mean.epi.eu.d.sub[q1.rois])
    eu.q2 <- mean(mean.epi.eu.d.sub[q2.rois])
    eu.q3 <- mean(mean.epi.eu.d.sub[q3.rois])
    eu.q4 <- mean(mean.epi.eu.d.sub[q4.rois])
    
    df.fc.sub <- data.frame(Quantile = c("Q1","Q2","Q3","Q4"),
                            Quantile.num = 1:4,
                            FC.quantile = c(fc.q1,fc.q2,fc.q3,fc.q4),
                            mean.sv2a = c(w.q1,w.q2,w.q3,w.q4),
                            EU.dist = c(eu.q1,eu.q2,eu.q3,eu.q4))
    
  })
  
  df.stat.sub.cluster <- merge(df.ad.epicenter.cluster,df.stat.sub.cluster, by ="SubID")
  mod.lmer <- lmer(data = df.stat.sub.cluster,mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB +EU.dist+(1|SubID)) ;print(summary(mod.lmer))
  
  # plot
  y.range <- range(df.stat.sub.cluster$mean.sv2a)
  p <- ggplot(data= df.stat.sub.cluster)+
    aes(x=Quantile, y = mean.sv2a)+
    geom_point( size=2, color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
    geom_line(color = "#DF8F44FF", alpha=0.5,aes(group=SubID))+
    geom_smooth(method = "lm",color = "grey55",aes(x=Quantile.num, y = mean.sv2a))+
    theme(panel.grid.major = element_blank(),
          panel.grid.minor = element_blank(),
          panel.background = element_blank(),
          panel.border = element_rect(color= " black", fill = NA),
          axis.text.x = element_text(size=10),
          axis.text.y = element_text(size=10))+
    xlab(str_wrap("Functional connectivity quantile to subject-specific epicenter", width = 40))+
    ylab(str_wrap("SV2A-PET W-score", width = 30))+
    ggtitle(paste0("Subject-level epicenter connectivity-based prediciton for subtype ",i))
  print(p)
  
}
## [1] "Subject-level analysis for cluster 1"
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +  
##     (1 | SubID)
##    Data: df.stat.sub.cluster
## 
## REML criterion at convergence: -134.3
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -2.50254 -0.71621 -0.06735  0.63409  2.83313 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.001155 0.03399 
##  Residual             0.024739 0.15729 
## Number of obs: 208, groups:  SubID, 52
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)      -0.522823   0.418079 200.667280  -1.251    0.213
## FC.quantile       0.040789   0.040758 155.784056   1.001    0.318
## SEXmale           0.026953   0.025238  47.039425   1.068    0.291
## Age.SDM8         -0.001613   0.001417  47.199081  -1.139    0.261
## DX.ABDementia_1   0.006228   0.027505  47.619268   0.226    0.822
## DX.ABMCI_1       -0.015792   0.035193  47.204933  -0.449    0.656
## EU.dist           0.085504   0.054673 194.805403   1.564    0.119
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.052                                   
## SEXmale      0.027 -0.008                            
## Age.SDM8    -0.207 -0.012 -0.056                     
## DX.ABDmnt_1  0.057  0.011 -0.241  0.023              
## DX.ABMCI_1  -0.042  0.003 -0.258  0.299  0.399       
## EU.dist     -0.965  0.056 -0.029 -0.053 -0.089 -0.054
## `geom_smooth()` using formula = 'y ~ x'
## [1] "Subject-level analysis for cluster 2"
## boundary (singular) fit: see help('isSingular')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +  
##     (1 | SubID)
##    Data: df.stat.sub.cluster
## 
## REML criterion at convergence: -25.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -2.7761 -0.5048 -0.1646  0.6790  2.3266 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.00000  0.000   
##  Residual             0.02788  0.167   
## Number of obs: 72, groups:  SubID, 18
## 
## Fixed effects:
##                  Estimate Std. Error        df t value Pr(>|t|)
## (Intercept)     -0.558297   0.812288 65.000000  -0.687    0.494
## FC.quantile     -0.034249   0.070572 65.000000  -0.485    0.629
## SEXmale          0.045913   0.044589 65.000000   1.030    0.307
## Age.SDM8         0.003379   0.003085 65.000000   1.095    0.277
## DX.ABDementia_1 -0.017758   0.050562 65.000000  -0.351    0.727
## DX.ABMCI_1       0.017808   0.052323 65.000000   0.340    0.735
## EU.dist          0.040302   0.104171 65.000000   0.387    0.700
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile  0.030                                   
## SEXmale     -0.108 -0.027                            
## Age.SDM8    -0.298 -0.006  0.265                     
## DX.ABDmnt_1  0.140 -0.015  0.233 -0.221              
## DX.ABMCI_1   0.127  0.011 -0.118  0.138  0.234       
## EU.dist     -0.961 -0.028  0.008  0.026 -0.113 -0.187
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## `geom_smooth()` using formula = 'y ~ x'
## [1] "Subject-level analysis for cluster 3"
## boundary (singular) fit: see help('isSingular')

## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile + SEX + Age.SDM8 + DX.AB + EU.dist +  
##     (1 | SubID)
##    Data: df.stat.sub.cluster
## 
## REML criterion at convergence: -104
## 
## Scaled residuals: 
##      Min       1Q   Median       3Q      Max 
## -3.01054 -0.65480  0.06666  0.55215  2.93787 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.00000  0.0000  
##  Residual             0.02605  0.1614  
## Number of obs: 172, groups:  SubID, 43
## 
## Fixed effects:
##                   Estimate Std. Error         df t value Pr(>|t|)
## (Intercept)     -6.747e-02  4.523e-01  1.650e+02  -0.149    0.882
## FC.quantile      2.070e-02  4.444e-02  1.650e+02   0.466    0.642
## SEXmale          2.830e-02  2.724e-02  1.650e+02   1.039    0.300
## Age.SDM8         4.788e-04  1.443e-03  1.650e+02   0.332    0.740
## DX.ABDementia_1 -1.251e-02  3.453e-02  1.650e+02  -0.362    0.717
## DX.ABMCI_1      -3.053e-02  3.113e-02  1.650e+02  -0.981    0.328
## EU.dist          2.993e-03  6.196e-02  1.650e+02   0.048    0.962
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt SEXmal A.SDM8 DX.ABD DX.ABM
## FC.quantile -0.098                                   
## SEXmale     -0.009 -0.010                            
## Age.SDM8    -0.024 -0.034  0.114                     
## DX.ABDmnt_1  0.000  0.002  0.326 -0.037              
## DX.ABMCI_1   0.032 -0.005 -0.035 -0.100  0.549       
## EU.dist     -0.969  0.104 -0.055 -0.214 -0.040 -0.048
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
## `geom_smooth()` using formula = 'y ~ x'

Analysis 5: epicenter connectivity-based prediciton of longitudinal change rates of SV2A wsocre

# define function 
get.epi.pred.beta.followup.shuffled <- function(epicenter.rois, mean.sv2a.w, mean.sv2a.w.cr){
  mapply(function(j){
    fc.mat.shuffled <- shuffled.fc[[j]]
    
    epi.fc.mat.shuffled <- as.matrix(fc.mat.shuffled[epicenter.rois,])
    epi.fc.mat.shuffled[which(epi.fc.mat.shuffled==Inf)] <- NA
    mean.epi.fc.mat.shuffled <- colMeans(epi.fc.mat.shuffled,na.rm = T)
    
    df.test.ad.shuffled <- data.frame(epi.fc = mean.epi.fc.mat.shuffled,
                                      sv2a.wscore=mean.sv2a.w,
                                      sv2a.wscore.cr = mean.sv2a.w.cr)
    mod.lm.shuffled <- lm(data = df.test.ad.shuffled,sv2a.wscore~epi.fc+mean.sv2a.w.cr)
    mod.lm.shuffled <- summary(lm.beta(mod.lm.shuffled))
    beta.val <- mod.lm.shuffled$coefficients[2,2]
    beta.val
  }, 1:1000)

}

# group-level analyses
df.ad.followup <- subset(df.ad.epicenter, !is.na(CR.W.SDM8.V1))
mean.w.bl <- colMeans(df.ad.followup[,paste0("W.SDM8.V",1:200)])
mean.w.cr <- colMeans(df.ad.followup[,paste0("CR.W.SDM8.V",1:200)])
epicenter.bl <- get.epicenter.by.rank(mean.w.bl)
epicenter.bl <- which(epicenter.bl==1)

epi.fc.mat.bl <- fc.mat[epicenter.bl,]
mean.epi.fc.bl <- colMeans(epi.fc.mat.bl,na.rm = T)

beta.shuffled.followup <- get.epi.pred.beta.followup.shuffled(epicenter.bl,mean.w.bl,mean.w.cr)

df.test.followup <- data.frame(epicenter.fc = mean.epi.fc.bl,
                               sv2a.wscore.cr = mean.w.cr,
                               sv2a.wscore.bl = mean.w.bl)

mod.lm <- lm(data = df.test.followup,sv2a.wscore.cr~epicenter.fc+sv2a.wscore.bl)
mod.lm <- summary(lm.beta(mod.lm));print(mod.lm)
## 
## Call:
## lm(formula = sv2a.wscore.cr ~ epicenter.fc + sv2a.wscore.bl, 
##     data = df.test.followup)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.279083 -0.066378  0.003909  0.073033  0.245720 
## 
## Coefficients:
##                Estimate Standardized Std. Error t value Pr(>|t|)    
## (Intercept)    -0.25229           NA    0.00753 -33.505  < 2e-16 ***
## epicenter.fc   -0.07310     -0.20377    0.02503  -2.920  0.00391 ** 
## sv2a.wscore.bl  0.01209      0.02965    0.02845   0.425  0.67139    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.1056 on 197 degrees of freedom
## Multiple R-squared:  0.04299,    Adjusted R-squared:  0.03327 
## F-statistic: 4.425 on 2 and 197 DF,  p-value: 0.01319
beta.val <- mod.lm$coefficients["epicenter.fc","Standardized"]
p.val.shuffled <- length(which(beta.shuffled.followup< beta.val))/length(beta.shuffled.followup)
# subject-level analyses
beta.sub.followup <- ddply(df.ad.followup,.(SubID), function(x){
  mean.w.sub.cr <- as.vector(as.matrix(x[,paste0("CR.W.SDM8.V",1:200)]))
  mean.w.sub.bl <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
  epicenters.bl.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
  
  epi.fc.mat.sub <- fc.mat[epicenters.bl.sub,]
  mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
  
  df.sub.tmp <- data.frame(mean.w.sub.cr = mean.w.sub.cr,
                           mean.w.sub.bl = mean.w.sub.bl,
                           mean.fc.sub = mean.epi.fc.sub)
  beta.sub <- coefficients(lm.beta(lm(data=df.sub.tmp,mean.w.sub.cr~mean.fc.sub+mean.w.sub.bl)))[2]
  c(beta.sub)
})

beta.sub.followup$variable <- "beta"
t.mod <- t.test(beta.sub.followup$mean.fc.sub)
t.mod
## 
##  One Sample t-test
## 
## data:  beta.sub.followup$mean.fc.sub
## t = -1.0008, df = 16, p-value = 0.3318
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.04683032  0.01679366
## sample estimates:
##   mean of x 
## -0.01501833

Analysis 6: effect of plasma ptau concentration on connectivity-mediated synaptic loss propagation

# determine subject-specific epicenters for all participants 
df.epicenters.subject.all <- ddply(df,.(SubID), function(x){
      mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
      epicenters.sub <- get.epicenter.by.rank(mean.w.sub)
      epicenters.sub
    })
colnames(df.epicenters.subject.all)[2:201] <- paste0("Epicenter.bin.V",1:200)
df.epicenter <- merge(df,df.epicenters.subject.all, by ="SubID")

# derive subject-level connectivity-based synaptic loss
beta.sub.all <- ddply(df.epicenter,.(SubID), function(x){
  mean.w.sub <- as.vector(as.matrix(x[,paste0("W.SDM8.V",1:200)]))
  epicenters.sub <- which(x[,paste0("Epicenter.bin.V",1:200)]==1)
  
  epi.fc.mat.sub <- fc.mat[epicenters.sub,]
  mean.epi.fc.sub <- colMeans(epi.fc.mat.sub,na.rm = T)
  
  df.sub.tmp <- data.frame(mean.w.sub = mean.w.sub,
                           mean.fc.sub = mean.epi.fc.sub)
  beta.sub <- coefficients(lm.beta(lm(data=df.sub.tmp,mean.w.sub~mean.fc.sub)))[2]
  c(beta.sub)
})
colnames(beta.sub.all)[2] <- "connectivity.SV2A.beta"
t.test(beta.sub.all$connectivity.SV2A.beta)
## 
##  One Sample t-test
## 
## data:  beta.sub.all$connectivity.SV2A.beta
## t = -0.1161, df = 149, p-value = 0.9077
## alternative hypothesis: true mean is not equal to 0
## 95 percent confidence interval:
##  -0.01150906  0.01023165
## sample estimates:
##     mean of x 
## -0.0006387068
df.all.ptau181 <- merge(df.epicenter, beta.sub.all, by ="SubID")
df.ad.ptau181  <- subset(df.all.ptau181, AV45.bin==1)


# linear regression
mod.all <- lm(data = df.all.ptau181,connectivity.SV2A.beta~ log(pTau181)+ Age.SDM8 +AV45.global.mean+ SEX + DX); summary(lm.beta(mod.all))
## 
## Call:
## lm(formula = connectivity.SV2A.beta ~ log(pTau181) + Age.SDM8 + 
##     AV45.global.mean + SEX + DX, data = df.all.ptau181)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.155639 -0.046858 -0.001261  0.043153  0.168463 
## 
## Coefficients:
##                    Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)       0.0247728           NA  0.0493242   0.502    0.616
## log(pTau181)      0.0048930    0.0719712  0.0057353   0.853    0.395
## Age.SDM8         -0.0002539   -0.0338735  0.0006301  -0.403    0.688
## AV45.global.mean  0.0020351    0.0215027  0.0080230   0.254    0.800
## SEXmale           0.0049042    0.0365165  0.0116550   0.421    0.675
## DXDementia       -0.0115818   -0.0750130  0.0135697  -0.854    0.395
## DXMCI            -0.0145663   -0.0888641  0.0144731  -1.006    0.316
## 
## Residual standard error: 0.06822 on 143 degrees of freedom
## Multiple R-squared:  0.01597,    Adjusted R-squared:  -0.02532 
## F-statistic: 0.3867 on 6 and 143 DF,  p-value: 0.8866
mod.ad <- lm(data = df.ad.ptau181,connectivity.SV2A.beta~ log(pTau181)+ Age.SDM8 +AV45.global.mean+ SEX + DX); summary(lm.beta(mod.ad))
## 
## Call:
## lm(formula = connectivity.SV2A.beta ~ log(pTau181) + Age.SDM8 + 
##     AV45.global.mean + SEX + DX, data = df.ad.ptau181)
## 
## Residuals:
##       Min        1Q    Median        3Q       Max 
## -0.137337 -0.045001 -0.002558  0.045269  0.165296 
## 
## Coefficients:
##                    Estimate Standardized Std. Error t value Pr(>|t|)
## (Intercept)      -0.0111801           NA  0.0570300  -0.196    0.845
## log(pTau181)      0.0079342    0.1150652  0.0066298   1.197    0.234
## Age.SDM8          0.0002119    0.0281121  0.0007275   0.291    0.771
## AV45.global.mean  0.0084790    0.0891837  0.0092927   0.912    0.364
## SEXmale           0.0084778    0.0634181  0.0134060   0.632    0.528
## DXDementia       -0.0175174   -0.1240537  0.0151090  -1.159    0.249
## DXMCI            -0.0184616   -0.1246824  0.0161495  -1.143    0.256
## 
## Residual standard error: 0.06756 on 106 degrees of freedom
## Multiple R-squared:  0.03807,    Adjusted R-squared:  -0.01637 
## F-statistic: 0.6993 on 6 and 106 DF,  p-value: 0.6508
# interaction effect between ptau181 and epicenter functional connectivity
## test for pooled sample
df.stat.plasma.ptau181.interaction <- subset(df.stat.sub, !is.na(pTau181))
df.stat.plasma.ptau181.interaction$pTau181.bin <- ifelse(df.stat.plasma.ptau181.interaction$pTau181>= median(df.stat.plasma.ptau181.interaction$pTau181, na.rm = T),1,0)

mod.lmer <- lmer(data = df.stat.plasma.ptau181.interaction,
                 mean.sv2a ~ FC.quantile*log(pTau181) + SEX + Age.SDM8 + DX + AV45.global.mean + EU.dist +(1|SubID))
## boundary (singular) fit: see help('isSingular')
summary(mod.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile * log(pTau181) + SEX + Age.SDM8 + DX +  
##     AV45.global.mean + EU.dist + (1 | SubID)
##    Data: df.stat.plasma.ptau181.interaction
## 
## REML criterion at convergence: -309.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0100 -0.6880  0.0005  0.5983  3.0817 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.00000  0.0000  
##  Residual             0.02587  0.1608  
## Number of obs: 452, groups:  SubID, 113
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)              -3.404e-01  2.831e-01  4.420e+02  -1.203    0.230  
## FC.quantile               4.635e-02  4.124e-02  4.420e+02   1.124    0.262  
## log(pTau181)              4.062e-03  7.900e-03  4.420e+02   0.514    0.607  
## SEXmale                   2.954e-02  1.597e-02  4.420e+02   1.850    0.065 .
## Age.SDM8                 -5.000e-04  8.713e-04  4.420e+02  -0.574    0.566  
## DXDementia               -1.037e-03  1.804e-02  4.420e+02  -0.057    0.954  
## DXMCI                    -1.898e-02  1.928e-02  4.420e+02  -0.985    0.325  
## AV45.global.mean         -1.011e-03  1.106e-02  4.420e+02  -0.091    0.927  
## EU.dist                   5.008e-02  3.773e-02  4.420e+02   1.327    0.185  
## FC.quantile:log(pTau181)  2.473e-02  2.982e-02  4.420e+02   0.829    0.407  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt l(T181 SEXmal A.SDM8 DXDmnt DXMCI  AV45.. EU.dst
## FC.quantile -0.039                                                        
## log(pTa181)  0.071  0.012                                                 
## SEXmale      0.002  0.000  0.131                                          
## Age.SDM8    -0.120 -0.011 -0.017  0.060                                   
## DXDementia   0.063 -0.002  0.000  0.035 -0.065                            
## DXMCI        0.035 -0.008 -0.028 -0.175  0.049  0.424                     
## AV45.glbl.m -0.036  0.002  0.044  0.152 -0.057 -0.035  0.075              
## EU.dist     -0.971  0.042 -0.043 -0.042 -0.110 -0.077 -0.074  0.000       
## FC.q:(T181) -0.002  0.741  0.003  0.007  0.004 -0.005 -0.008  0.002  0.002
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
y.range <- range(df.stat.plasma.ptau181.interaction$mean.sv2a)

ggplot(data=df.stat.plasma.ptau181.interaction)+
  aes(x=Quantile ,y=mean.sv2a, color=factor(pTau181.bin))+
  geom_point(shape=19,alpha=0.5)+
  geom_line(alpha=0.15,aes(group=SubID))+
  geom_smooth(method = "lm",aes(x=Quantile.num , y = mean.sv2a, color =factor(pTau181.bin)), linewidth=0.5)+
  labs(x = "",y = "SV2A-PET W-score") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color= " black", fill = NA),
        axis.text.x = element_text(size=12))+
  scale_x_discrete()+
  scale_color_nejm(name="Plasma p-tau181",labels=c("<median",">median"))+
  guides(fill = "none")
## `geom_smooth()` using formula = 'y ~ x'

## test for abeta positive sample
df.stat.plasma.ptau181.interaction.ad <- subset(df.stat.sub, !is.na(pTau181)& AV45.bin==1)
df.stat.plasma.ptau181.interaction.ad$pTau181.bin <- ifelse(df.stat.plasma.ptau181.interaction.ad$pTau181>= median(df.stat.plasma.ptau181.interaction.ad$pTau181, na.rm = T),1,0)

mod.lmer <- lmer(data = df.stat.plasma.ptau181.interaction.ad,
                 mean.sv2a ~ FC.quantile*log(pTau181) + SEX + Age.SDM8 + DX +AV45.global.mean + EU.dist +(1|SubID))
## boundary (singular) fit: see help('isSingular')
summary(mod.lmer)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: mean.sv2a ~ FC.quantile * log(pTau181) + SEX + Age.SDM8 + DX +  
##     AV45.global.mean + EU.dist + (1 | SubID)
##    Data: df.stat.plasma.ptau181.interaction.ad
## 
## REML criterion at convergence: -309.8
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.0100 -0.6880  0.0005  0.5983  3.0817 
## 
## Random effects:
##  Groups   Name        Variance Std.Dev.
##  SubID    (Intercept) 0.00000  0.0000  
##  Residual             0.02587  0.1608  
## Number of obs: 452, groups:  SubID, 113
## 
## Fixed effects:
##                            Estimate Std. Error         df t value Pr(>|t|)  
## (Intercept)              -3.404e-01  2.831e-01  4.420e+02  -1.203    0.230  
## FC.quantile               4.635e-02  4.124e-02  4.420e+02   1.124    0.262  
## log(pTau181)              4.062e-03  7.900e-03  4.420e+02   0.514    0.607  
## SEXmale                   2.954e-02  1.597e-02  4.420e+02   1.850    0.065 .
## Age.SDM8                 -5.000e-04  8.713e-04  4.420e+02  -0.574    0.566  
## DXDementia               -1.037e-03  1.804e-02  4.420e+02  -0.057    0.954  
## DXMCI                    -1.898e-02  1.928e-02  4.420e+02  -0.985    0.325  
## AV45.global.mean         -1.011e-03  1.106e-02  4.420e+02  -0.091    0.927  
## EU.dist                   5.008e-02  3.773e-02  4.420e+02   1.327    0.185  
## FC.quantile:log(pTau181)  2.473e-02  2.982e-02  4.420e+02   0.829    0.407  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) FC.qnt l(T181 SEXmal A.SDM8 DXDmnt DXMCI  AV45.. EU.dst
## FC.quantile -0.039                                                        
## log(pTa181)  0.071  0.012                                                 
## SEXmale      0.002  0.000  0.131                                          
## Age.SDM8    -0.120 -0.011 -0.017  0.060                                   
## DXDementia   0.063 -0.002  0.000  0.035 -0.065                            
## DXMCI        0.035 -0.008 -0.028 -0.175  0.049  0.424                     
## AV45.glbl.m -0.036  0.002  0.044  0.152 -0.057 -0.035  0.075              
## EU.dist     -0.971  0.042 -0.043 -0.042 -0.110 -0.077 -0.074  0.000       
## FC.q:(T181) -0.002  0.741  0.003  0.007  0.004 -0.005 -0.008  0.002  0.002
## optimizer (nloptwrap) convergence code: 0 (OK)
## boundary (singular) fit: see help('isSingular')
y.range <- range(df.stat.plasma.ptau181.interaction.ad$mean.sv2a)

ggplot(data=df.stat.plasma.ptau181.interaction.ad)+
  aes(x=Quantile ,y=mean.sv2a, color=factor(pTau181.bin))+
  geom_point(shape=19,alpha=0.5)+
  geom_line(alpha=0.15,aes(group=SubID))+
  geom_smooth(method = "lm",aes(x=Quantile.num , y = mean.sv2a, color =factor(pTau181.bin)), linewidth=0.5)+
  labs(x = "",y = "SV2A-PET W-score") +
  theme(panel.grid.major = element_blank(),
        panel.grid.minor = element_blank(),
        panel.background = element_blank(),
        panel.border = element_rect(color= " black", fill = NA),
        axis.text.x = element_text(size=12))+
  scale_x_discrete()+
  scale_color_nejm(name="Plasma p-tau181",labels=c("<median",">median"))+
  guides(fill = "none")
## `geom_smooth()` using formula = 'y ~ x'